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Abstract 

We investigate the thermal equilibrium properties of kinks in a classical <1>^ field theory in 
1 + 1 dimensions. From large scale Langevin simulations we identify the temperature below 
which a dilute gas description of kinks is valid. The standard dilute gas/WKB description 
is shown to be remarkably accurate below this temperature. At higher, "intermediate" tem- 
peratures, where kinks still exist, this description breaks down. By introducing a double 
Gaussian variational ansatz for the eigenfunctions of the statistical transfer operator for 
the system, we are able to study this region analytically. In particular, our predictions for 
the number of kinks and the correlation length are in agreement with the simulations. The 
double Gaussian prediction for the characteristic temperature at which the kink description 
ultimately breaks down is also in accord with the simulations. We also analytically calcu- 
late the internal energy and demonstrate that the peak in the specific heat near the kink 
characteristic temperature is indeed due to kinks. In the neighborhood of this temperature 
there appears to be an intricate energy sharing mechanism operating between nonlinear 
phonons and kinks. 
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I. Introduction 

The equilibrium and nonequilibrium statistical mechanics of solitons, solitary 
waves, and other coherent structures in nonlinear systems has been a subject of study 
for some time Recent interest has been fueled by new applications not only in 
condensed matter physics but also by potential applications in particle physics 

(sphalerons) and cosmology (domain walls, baryogenesis) [Q. 

In this paper we focus on the classical equilibrium statistical mechanics of solitary 
wave ( "kink" ) solutions of a tachyonic mass, 0^ ( "double well" ) field theory in 1 + 1 
spacetime dimensions with the Lagrangian density 

L = l {dt^f - \ {d^^f + \mH' - \^^^ (1) 

and the corresponding equation of motion 

= dl^^ + - A$l (2) 

This model is of direct relevance to the study of displacive phase transitions and 
magnetic spin chains 0. Moreover, the behavior of this model is representative of a 
large class of soliton-bearing systems. Finally, it has the added advantage of being 
amenable to both theoretical analysis and numerical simulation. 

The kink solution of the equation of motion (0) is the localized field configuration 
that interpolates between the two asymptotic field potential energy minima at $ = 
±$0 = ^fn/ -\/A. Kink solutions of the field equations are inaccessible to perturbation 
theory: A statistical mechanics of kinks is nevertheless still possible, partly because 
they are localized objects. The static kink solution centered aX x = Xq is 

$fc = -^tanh I -^(x - Xo) I , (3) 

and the negative of the kink solution is the anti-kink. Time dependent solutions are 
easily found by Lorentz boosting the static solution. The energy density of the static 
kink (or anti-kink) is 

e,(x) = -sech (4) 
from which the total energy of an isolated kink, or the kink mass, is obtained: 
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In numerical simulations it is customary to employ a dimensionless form of the 
theory given by the transformations: 
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where = m? / K. Applying these transformations, we have 

H = ma^H 



where H is the original field Hamiltonian, and 
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The new field satisfies the equation of motion 
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We enforce = PH by introducing a new scaled temperature P = [3 /{ma?). 

The statistical mechanics of kinks in this system has been studied by two ap- 
proaches. In the first, and phenomenological, approach one assumes that the kinks 
and the field fluctuations about the asymptotic field minima ("phonons") may be 
treated as weakly interacting elementary excitations. Provided that the kink density 
is low (the dilute gas approximation), the canonical partition function can be found 
by standard methods Alternatively, as shown by Krumhansl and Schrieffer 

(KS) 1^, building on earlier work of Sears, Scalapino, and Ferrell |]10[, it is possible to 
calculate the partition function, in principle exactly, by exploiting a transfer operator 
technique. KS showed that in the low temperature (dilute gas) limit the partition 
function naturally factorizes into two contributions both having counterparts in the 
phenomenological theory; a tunneling term which they were able to identify with the 
kink contribution, and the remainder which they identified as the linearized phonons. 
The approach of KS was further refined and extended to a wider class of systems 
by Currie, Krumhansl, Bishop, and Trullinger In this work, interactions of kinks 
with linearized phonons were taken into account, leading to substantial corrections 
to the results of KS. 

The key result of these efforts is the prediction that, below a certain temperature, 
the spatial density of kinks is 



Uk oc ^Ek(3exp{-Ek(3). 
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A related prediction is that the field correlation length A defined by 

(<l>(0)<l>(x)) ~ exp(- \x\ /A) (12) 
exhibits an exponential temperature dependence p, 

at low temperatures. 

Computer simulations to verify these results date back to Ref. |TT[ where only a 
quahtative agreement was found. Recent work |]12||13|[1^ has led to more detailed 
comparisons, however significant discrepancies have been reported. This has led to 
theoretical speculation 115||16| regarding possible corrections to the dilute gas theory 
of kinks. It has been suggested that these discrepancies are due to finite size effects 
and phonon dressing of the bare kink energy (breather contributions to the free energy 
may also be significant |jl8|). As was discussed in Ref. |jl9[, the earlier simulations 
were not carried out at low enough temperatures: Nevertheless the authors interpreted 
their results in terms of WKB formulae that are simply not a valid description over 
the range of temperatures they had studied. In this intermediate temperature regime, 
there is no unique characterization of what constitutes a kink and how to differentiate 
it from a nonlinear phonon. By going to low enough temperatures where the dilute 
gas results are valid, we are also able to rule out an earlier claim of substantial phonon 
dressing even at these temperatures fill . 

We have numerically studied the equilibrium statistical mechanics of kinks in the 

model by implementing a Langevin code on a massively parallel computer. To un- 
derstand our results in the high and intermediate temperature region not susceptible 
to a dilute gas analysis, we have used a nonperturbative double Gaussian wave func- 
tion approximation in the quantum mechanical problem for determining the eigen- 
values of the transfer operator. The study of the intermediate temperature regime 
is important because kink contributions to some thermodynamic quantities (e.g., the 
specific heat) of the system may be dominant precisely in this region. Numerical 
calculations of the partition function, while certainly valuable, are by themselves not 
sufficient. Such calculations, for example, cannot explain the nature of the peak in 
the specific heat nor can they enable one to extract which effects are due to kinks 
and which effects due to phonons. 

Stated in brief, our results are: (1) the dilute gas predictions for the kink density 
and the correlation length are very accurate below a certain (theoretically estimable) 
temperature, (2) above this temperature the Gaussian results for the kink number 
and correlation length agree with the simulations, (3) kinks are found to "disap- 
pear" into the thermal phonon background above a characteristic temperature, in 
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good agreement with our theoretical prediction, (4) our Gaussian approximation ac- 
curately describes the classical single point field distribution function at high and 
intermediate temperatures where the dilute gas (WKB) approximation breaks down, 
(5) the internal energy and the specific heat calculated in double Gaussian approxima- 
tion show an interesting energy sharing process between kinks and nonlinear phonons 
in an intermediate temperature range below the characteristic temperature at which 
kinks appear: a peak in the specific heat in this temperature range is shown to be 
due essentially to kinks. 

The rest of this paper is organized as follows. We describe our computer simulation 
and numerical techniques in Sec. II. Sec. Ill is a brief review of the transfer operator 
formalism and the standard WKB results for the partition function. Here we also 
introduce a new, "effective" definition of the kink number that takes into account 
an averaging over the thermal phonon length scale. In Sec. IV we introduce the 
double Gaussian variational method and use it to find approximate eigenfunctions and 
eigenvalues of the transfer operator and also to compute thermodynamic quantities 
such as the specific heat. Finally, we conclude in Sec. V with a discussion of our 
results and of directions for future work. 
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II. The Simulations 

The Langevin equation for the dimensionless theory is 

dU = dl,<P - vdt<P + 0(1 - 0') + Fix, i). (14) 

To guarantee an approach to equilibrium, the Gaussian, white noise F and the vis- 
cosity rj are related via the fluctuation-dissipation theorem: 

{F{x,t)F{y,s)) =2r]p-H{x-mt-s). (15) 

Since in this paper we are interested only in the equilibrium properties of the system 
which are independent of the viscosity, we flxed t] = 1 for all the simulations. (It was 
verifled that the results did not depend on the value of t].) 

We carried out numerical simulations of this Langevin equation using a standard 
flrst order Euler differencing technique. The second order in time Langevin equation 



(p!^) was written as the two flrst order flnite difference equations 

v{i+e) = v{i) + e[(^^{x + 6) + ^{x-6)-2^{x)) 

-r]v{i) + 0(x)(l - 02(x)) + F{x, t) 
0(t + e) = ^{t) + ev{t + e) (16) 

where the "velocity," v = dt4>. The flnite differencing was implemented with a time 
step e = 0.02, and a lattice spacing 6 = 0.5. The spacetime noise was generated by 
summing random variables uniformly distributed over (—.5, .5). From the central 
limit theorem, in the limit that N —>■ (yo this sum should approach a Gaussian random 
variable with mean and variance \J N/Vl. For convenience we chose = 12: 
since this gave a noise normalized to a strength of unity, for each given temperature 
the normalization was straightforwardly determined by the fluctuation-dissipation 
relation. The results presented here were obtained on lattices consisting of 16384 
sites and we have checked their consistency with results from lattices of different 
sizes. The lattice volume is large enough that there are no discernible flnite size 
effects. 

Our system size is one to two orders of magnitude larger than that in most previous 
simulations. Large system sizes are necessary to get acceptable statistics at low 
temperatures. In many of the recent numerical simulations of this model, the system 
size was such that only a few kinks would appear at the lowest temperatures studied. 
For temperatures in the range where the WKB theory is valid, fewer than one kink 
would appear on systems of the sizes used in Refs. [|12[ |jl3|] |14 . 



The lattice systems were evolved from random initial conflgurations to equilib- 
rium. The length of time necessary to ensure equilibrium increased with inverse 
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temperature. For (3 = 8 the time required was approximately 10'' time steps, and for 
the highest temperatures, less than 10^ steps. 

Two quantities of interest reported here are the kink number and the field corre- 
lation length. To compute the kink number, we need an operational way to identify 
kinks, even though there is an exact kink solution available theoretically. We therefore 
examined several possible definitions, all of which rely on a knowledge of the canonical 
kink size. From the classical solution for a kink centered at Xq, (j) = tanh[[x—XQ) / 
the kink scale is approximately 8 lattice units. Raw kink configurations are shown 
in Fig. I. At low enough temperatures (/3 > 5), kinks may be identified easily, however 
at higher temperatures this is clearly not the case. 

The naive approach is to simply count the number of zero-crossings of the field, 
since one may argue that these are the "tunneling events" which correspond to 
kinks. However, at higher temperatures there are zero-crossings due to thermal noise 
(phonons), and counting all zero- crossings would lead to a gross overestimation of the 
number of kinks. At high temperatures it is not possible to distinguish unambigu- 
ously between kinks and nonlinear phonons. However, the problem of the possible 
overcounting arises already at intermediate temperatures, where kinks are distinct. 
A possible solution is to use a smoothed field by either averaging or "block-spinning" 
the actual field configuration over a length of the order of the kink length scale. The 
latter approach was taken in previous simulations |T^0][|l4l. This solution is not 
without flaws either, as rapid fluctuations can still appear as kinks. We prefer to 
count kinks in the following way: at a particular time we first find all zero-crossings. 
To test the legitimacy of a given zero-crossing we check for zero-crossings one kink 
scale (8 lattice units) to its right and to its left. If no zero-crossings are found, we 
count it kink, otherwise not. 

The number of kinks is plotted against (3 in Fig. II. Above /5 ~ 6, the averaged 
field method and our method for counting kinks agree. Moreover, in this (low tem- 
perature) range, the dilute gas expression for the kink number (|TTD is in excellent 
agreement with the data. At elevated temperatures, there is a clear disagreement be- 
tween the two methods of counting kinks. The average field technique has the number 
of kinks monotonically increasing with temperature; whereas, in accord with intuition 
and the behavior of -P[0], the second technique clearly shows a reduction in the kink 
number at higher temperatures. (See Sec. IV for a detailed discussion.) Moreover, in 
this temperature regime the number of kinks computed with the smoothing method 
depends strongly on the smoothing scale. We conclude that for /3 < 6, an unambigu- 
ous number of kinks cannot be extracted with any confidence from the smoothing 
method. Unfortunately, this is precisely the temperature regime explored in previous 
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simulations. We will discuss an analytic definition of the kink number N^. in Sec. Ill 
and show that this quantity can be independent of the smoothing scale only at low 
temperatures. 

The correlation length was extracted from the field configurations by taking the 
inverse Fourier transform of the power spectrum and then fitting an exponential decay 
to the correlation function. The correlation length A is plotted against P in Fig. III. 
The WKB prediction is seen to hold for /3 > 6 while at higher temperatures (/? > 4), 
the double Gaussian approximation is in excellent agreement with the data. 

Our numerical results can be understood by using different analytical approaches 
in different temperature regimes. At low temperatures {f3 > 6) the WKB method 
turns out to be very accurate. At higher temperatures a double Gaussian variational 
method provides a good description of our data. Both analytical methods are con- 
veniently discussed in the framework of the transfer operator formalism, which we 
describe in the following section. 
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III. The Transfer Operator 

The key idea behind the transfer operator method is to transform the problem of 
finding the canonical partition function for a system to the exactly equivalent prob- 
lem of finding the eigenvalues of a certain integral operator. With some smoothness 
assumptions for the eigenfunctions of this transfer operator, it is then possible to 
show that the problem reduces to that of finding the energy eigenvalues of a related 
quantum mechanical problem. The method easily generalizes to the problem of cal- 
culating correlation functions. The presentation given here closely follows that of Ref. 

The canonical partition function for the Lagrangian ([^) is given by the functional 
integral 

Z = J D^Dtt exp(-/5iJ[$,7r]) (17) 

where tt is the canonical momentum of the field and H is the field Hamiltonian. The 
integral over the momenta is a trivial Gaussian integral. Writing 

Z = Zt^Z^, (18) 

we have 

obtained by discretizing the spatial lattice with lattice spacing A. Periodic boundary 
conditions are assumed and the total lattice size L = NA. We will be primarily 
interested in obtaining the configurational partition function, Z^. Skeletonizing the 
formal path integral for Z^ by first forming 



n 2 



A 



(20) 



we have. 



TV 

Z^ = ll d^, exp (-/3A/($,+i, $,)) (21) 

i=l 

where as just stated we have imposed a periodic boundary condition ($i = $7v+i), 
and 

= (22) 

The normalization ( P^ is fixed by requiring equipartition to hold for a free scalar 
field theory. 

Introducing the transfer operator specified by the eigenvalue problem 

J exp (-/3A/(<l>,+i, $,)) ^n(<f .) = exp (-/3AeJ ^n($.+i) (23) 



where e„ is the eigenvalue corresponding to the eigenfunction Substituting in 
(^), it follows from the completeness of the eigenfunctions that 

Z$ = ^e~^^^". (24) 

n 

In the thermodynamic limit, L — oo, Z$ is determined by the smallest eigenvalue of 
the transfer operator. In this limit, the free energy density is simply Fq, = cq. 

In an entirely analogous manner, it is possible to show that the correlation function 

($(x)$(0)) = ^ ( e-^^[*]$(x)$(0) 

= 5]|(*^„|$|^o)|'exp(-/3|a;|Ao„), (25) 

n 

where Ao„ = e„ — eo- At large values of x, the lowest excited state controls the 
behavior of the correlation function. 

The problem is now to compute the eigenvalues and eigenfunctions of the transfer 
operator. Assuming that the \l/„($i) are smooth functions of the $i, we Taylor expand 
^„($j) on the left hand side of ( pSj) in terms of the \l/„($j+i) on the right hand side. 
The integral on the left hand side can now be evaluated term by term. To leading 
order in A//3, ( p3|) becomes 



exp {-(3^Hq) = exp (-/3Ae„) ^„ (26) 
which is equivalent to the time independent Schrodinger equation 

HQ^n = en^n. (27) 

Dropping the indices on the fields, 

Hq = --m"^^ + -\^^ + Vo (28) 

is the quantum Hamiltonian for a particle in a double well potential with a temper- 
ature dependent energy shift 

Vq = In ( ^1 (29) 

arising from the normalization (^2]). The role of h is now played by (3^^, so that as 
the inverse temperature is made larger, WKB becomes a better approximation. For 
the dimensionless form of the theory (0), the quantum Hamiltonian (^) becomes 
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where Vq is just Vq with (3 replaced by /3. This is the form of the theory that we will 
use in the remainder of the paper. 

It is important to note that since ( p7[ ) was derived under certain smoothness 
requirements, the eigenf unctions of Hq or of Hq will not always be eigenf unctions of 
the transfer operator. This is the case at high temperatures where, for fixed A, the 
ratio A/j3 may be much larger than unity. However, for the intermediate and low 
temperature regimes which are of relevance to kink thermodynamics, ( p7D provides 
an adequate description. (If needed, the high temperature behavior can be studied 
via perturbation theory PU[.) 

An alternative way to view (|30|) is to introduce the scaled field = (3(f), in which 
case. 

Here we have omitted the contribution of Vq which, since it is just a shift, can be 
treated separately. 

At low temperatures the two wells are widely separated and the ground state 
energy is given by the oscillator ground state energy for one of the wells minus the 
tunnel-splitting term, usually calculated by the WKB method. One can then write 

(32) 

or, in terms of the free energies, 

F = Fosc + Ftunn + Vq (33) 

where, from a crude WKB estimate, 
and 

L 

^ tunn ^ 



^-e-^'^^ (35) 



where, the kink energy for the dimensionless form of the theory, = y8/9. The 
key result of KS is the realization that the tunneling contribution is associated with 
kinks and the sum of all the other contributions with phonons 

More generally one expects that there should be three qualitatively distinct tem- 
perature ranges. At very low temperatures the tunnel- splitting between two "classi- 
cally degenerate" states is very small and the WKB picture of kinks is valid. At higher 
temperatures the wave function overlap in the region between the wells is large enough 
that the WKB result for the energy difference between the ground and first excited 
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state is incorrect. One expects therefore that for temperatures in this "intermediate 
overlap" regime, the WKB approximation breaks down, but that kinks still exist as 
fairly well identifiable objects. Moreover, the tunnel-splitting between the two ground 
states can become comparable to the splitting between a ground state and the first 
excited state inside each well (this, however, has no effect on the correlation length, 
and as seen in Fig II, makes only a minor correction to the kink number). Finally, at 
still higher temperatures when the ground state energy in each well becomes larger 
than the classical barrier height, the kink description fails since the overlap is of order 
unity and the kink and phonon energy scales become indistinguishable. These three 
temperature regimes are in fact clearly seen in our numerical data (Fig. I) as well as 
in the analytic results from the double Gaussian approximation of Sec. IV (Fig. V). 

Let us first describe the expectation values for the correlator and the kink density 
based on the WKB approximation, which one expects to be valid in the low temper- 
ature regimeQ. An accurate calculation in this limit (see, e.g., Ref. |§]) yields, for the 
energy difference between the ground and first excited state, 

/op \ 1/2 

^oi=4f-^j exp{-Ej). (36) 

The WKB result for the correlation function (using only the first two states) is 

{(j){0)(j){x)) ~ e-l^l/^ (37) 



where the correlation length 



A = (3^ 



1 [W 
4 



exp 



{EkP) . (39) 



This is plotted in Fig. Ill and is seen to fit the numerical data very well for (3 > 6. 

Another quantity of interest is the number of kinks. In fact there appears to be no 
unambiguous way to define this quantity. The number of kinks is usually calculated 
in the phenomenological approach by working in the grand canonical ensemble (see, 
e.g., Ref. IP). We now introduce a new way of defining an effective kink number by 
working only with the original field variable. We begin by identifying three length 



^From now on we perform all calculations in the thermodynamic limit L ^ oo. The large size of 
the system in our simulation (16384 sites) rules out finite size effects: from (p^), it is clear that the 
thermodynamic limit is correctly found in the simulation provided that LAoi/3 3> 1, or in terms of 
the correlation length, when L ^ X. At the lowest temperature studied here L/X > 40. Since this is 
the worst case result, the condition L X is very well satisfied over the full range of the simulation. 
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scales in the problem: the averaging scale, a length scale long enough such that 
the averaged field 

0a(x) = r / " dy4>{x + y) (40) 

ta -I -la/2 

is smooth on scales of the order of the kink length. Calling the typical kink separation, 
/fc, and with such that < << Z^, we define 

\2 \ 




0,2) - {UO)MLs))] (41) 



to be the kink number density. This definition comes from just counting the number of 
zero-crossings of the smoothed field under an assumption that the typical separation 
between the kinks is much larger than and Lg. The normalization factor ±0o is the 
asymptotic field value away from a kink. The angular brackets denote a sampling over 
the whole lattice (divided into blocks of length Lg). The factor 1/2 compensates for 
the fact that anti-kinks are also counted (the number of kinks is equal to the number 
of anti- kinks). Note that our definition of kink number is sensible only when there is 
an appropriate separation of length scales, e.g., it is not valid at high temperatures 
when the kink separation is very small (where even the notion of a kink is not well 
defined). Even at intermediate temperatures, A^^ depends on the smearing length. 
Such a dependence has indeed been noted in numerical simulations . 

To understand the errors in estimating the number of kinks from (^T]) we first 
consider the case of fixed and varying Lg. If, on the one hand, Lg is chosen too 
large {Lg > Ik) then there is an undercounting of the total number of kinks since 
there will be many instances of having more than one kink in a block of length Lg. 
If, on the other hand, Lg is too small (i.e., smaller than the kink size) then there 
will again be an undercounting since in any given block of size Lg, |0a(O) — (t>a{Lg)\ < 
200- Therefore, for a given /„, one should maximize (|41|) with respect to Lg. Now 
considering la as variable, it is clear that for small la {la less than the kink size), (41) 
can be nonzero as long as nonnegligible phonon fiuctuations are present independent 
of whether kinks exist or not. At low temperatures one expects A^^ to be independent 
of la (see the discussion below) but in the intermediate temperature regime there 
should be a dependence on the averaging scale. 



An explicit formula for A^^ can be obtained by substituting (^) in (|25|). Then, 
for Lg > la, it is easy to show that 

m = ,5-^ E ^^^^ [pAoJa - 1 + e-^'^'^o" + e-^^=^- (l - cosHPAoJa))' • 

(42) 
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In the dilute gas regime we keep only the contribution of the first excited state in 



(H), use the WKB result (gl) for Aqi, and take the limit /? oo in (H). The total 
number of kinks is, for Lg > la, 

which, modulo the constant prefactor, is in agreement with the result of Ref. in 
this limit. An interesting point is that we have an explicit expression for the prefactor 
and therefore can test our formula directly against the simulations. 

Our predictions are compared with results from the simulations in Fig. II. (Note 



that the kink energy in the dimensionless form of the theory = y8/9 and that 
EkjS = Ekj3.) In the low temperature limit, one expects the kink number to be 
independent of Lg, and this is indeed true provided that ^ Ls (see (^3D). In this 
case, while the analytic result overestimates the number of kinks by a multiplicative 
factor of ~ 1.6, the functional form of the temperature dependence of the kink number 
is in very good agreement with the simulations. If one sets Ls = la in (0) then 
compared with the previous case, the kink number is reduced by a factor of 2/3, and 
is in complete agreement with the simulations. At present we do not understand 
the sensitivity of the theoretical formula to the ratio of la/ Lg in the low temperature 
regime. A discussion of the kink number density in the intermediate temperature 
regime will be postponed to Sec. IV, after we have described the double Gaussian 
approximation. 

The agreement of the WKB/dilute gas results with the simulations imply that 
there is no significant renormalization of the kink energy (beyond that due to lin- 
earized phonons) at low temperatures (i.e., /3 > 6). This is in disagreement with the 
simulations carried out in Ref. [|r^ but in agreement with the theory of Ref. 0. 
Whether there is or not such a renormalization of the kink mass at intermediate tem- 
peratures is difficult to analyze as at these temperatures the effects due to nonlinear 
phonons and kinks are hard to disentangle. A good example of this is the behavior 
of the kink number versus (3 (Fig. II). Both numerical and analytic results show that 
the ambiguity in the very notion of a kink number density is such as to rule out any 
estimation of the kink mass from the data at intermediate temperatures. 

In the regime where WKB fails, one can compare the simulations of the kink sys- 
tem with numerical solutions for the energy eigenvalues of the Hamiltonian Hq. How- 
ever, one would like to have a simple analytical method for predicting the measured 
quantities. To this end we implement the double Gaussian variational approxima- 
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tion in the effective quantum mechanical problem. This approximation is an order of 
magnitude more accurate than the simple Gaussian approximation for this problem 
and correctly accounts for the reduction of energy due to the overlap between the 
wave functions in the two wells. It also allows a natural decomposition of phonon 
and kink contributions in terms of separate "diagonal" and "overlap" contributions 
to the ground state energy. The double Gaussian approximation should be reliable 
at intermediate and high temperatures, where the overlap between the two ground 
states is substantial (at low temperatures, when the overlap is very sensitive to the 
form of the "tail" of the wave function, WKB is more accurate). 
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IV. The Double Gaussian Approximation 

The Gaussian approximation is a well-known nonperturbative variational method 
for calculating the ground state energy and effective potentials in quantum mechan- 
ics and quantum field theory As discussed in the previous section, the transfer 



operator technique reduces the classical statistical mechanics of a field theory in 1 + 1 
dimensions to ordinary quantum mechanics. To apply the standard Gaussian approx- 
imation is simple: one assumes that the ground state wave function is a Gaussian 
with width ^l'^^"^, and centered at the point (pQ, 




^g{<Po, fi) = \/ r ( - <^o)' ) • (44) 
Next one computes the energy 

Vg{^o,^) = {^g\Hq\^g), (45) 

and minimizes with respect to Q to yield the Gaussian effective potential Vg(0o)- The 
global minimum of Vg'(0o) is the ground state energy as calculated in the Gaussian 
approximation. 

In our case, with the Hamiltonian given by (0), it is easy to show that 



Minimization of ( ^61) with respect to yields the cubic equation 




3 3 



0, (47) 



the largest positive root of which is of interest. Substitution of (^7]) in (|46|) enables 
us to write 

VcW = -^^„^ + jii^„* + in-^i,. (48) 

Note that 0o = is always a local minimum of Vg(0o)- Above a certain temperature, 
it becomes the global minimum and one has a behavior reminiscent of a first order 
phase transition in the effective potential (Fig. IV). Numerically, this transition 
occurs at /3 = 3.12 and one might interpret this as predicting a disappearance of 
kinks at this temperature. This is not quite correct, however, as we will soon see. 

A problem with the Gaussian approximation is that it does not account for tun- 
neling. We know intuitively that for the double well Hamiltonian, at large well sepa- 
ration, the ground state should be described by a superposition of two localized wave 
packets, one in each well. These wave packets may be taken to be Gaussians. With 
such a wave function there will be cross-terms in (^^1 Hq |^) arising from overlaps 
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between the individual Gaussians. Correctly accounting for these overlap terms is 
clearly essential for our problem. 

We therefore modify the trial wave function by taking it as a superposition of two 
Gaussians: 



DG 



N 



where the normalization factor 



1 + 



(49) 



(50) 



It is then a simple matter to show that 



V^DG(0O,fi) 
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16p^ ^ I, 



(51) 



V 2 Al3^n 

The first term represents the "diagonal" contribution to the energy while the rest 
are cross-terms from the overlap or "off-diagonal" contribution. In the framework 
of the phonon-kink model these terms have a natural interpretation: the "diagonal" 
term corresponds to a contribution to the free energy due to phonons (including 
anharmonic corrections) while the overlap term is the kink contribution. In the low 
temperature limit this reproduces the WKB splitting of the free energy into phonons 
and kinks. 

In principle, one can minimize Vdg by following the same procedure as for Vg- 
Unfortunately, this yields a transcendental equation for Q which cannot be handled 
analytically. A way out is to first minimize (ETI) with respect to Q ignoring the 
overlap terms completely: this amounts to keeping Q as given by { ^7\) in (51). The 
minimization with respect to (pQ is then carried out by differentiating (^T]) with respect 



to (f)o and setting the result to zero. While this procedure appears to be rather a 
drastic simplification, it is still a major improvement over Vg both qualitatively and 
quantitatively. First, Vdg(0o) always has a local maximum at 0o = so that there is 
no abrupt change in the behavior of Vdg(0o) with temperature, consistent with the 
fact that there is no finite temperature phase transition in this model {Vdg and Vg 
are compared in Fig. IV). Second, it is an order of magnitude more accurate than Vg 
in estimating the ground state energy (worst case error ~ 1% compared to ~ 15%). 
Third, and very importantly, we can now estimate the energy of the higher excited 
states. 

To estimate the energy of the first excited state we simply minimize vj^Q = 
(\I'i| Hq |\E'i) where \E'i is the antisymmetric partner of "^dg'- 

1 

N 
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with A^^ = 2[1 — exp(— f2i</)]^^)]. It is to be stressed that the minimization of the energy 
with respect to has to be done independently from that with ^ dg- We find that 
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The difference between the absolute minima of vj^Q and Vdg is the tunnel-sphtting 
energy Aqi calculated in the double Gaussian approximation. We now investigate the 
low temperature limit of the double Gaussian approximation and contrast it with the 
WKB result. First, from (|47|), it follows immediately that at low temperatures, 



(54) 



Also, at low temperatures, 
(ETi we find, 



p. Substituting these last two results in ( pT| ) and 
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Polynomial corrections in powers of to the leading order result are due to anhar- 
monic corrections while the exponential terms represent the tunneling contribution. 
On subtracting ( |55D from (|56| ) the tunnel-splitting term follows immediately: 
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(57) 



which may be contrasted with the WKB result (^). In this regime, the WKB result 
is more accurate as the double Gaussian method tends to underestimate the overlap 

exp(— -\/2/^) in double Gaussian versus Aqi 



contribution (Aqi ~ exp(— v2/5) in double Gaussian versus Aqi ~ exp(— y8/9/5) 
in WKB). This is to be expected since in this regime the overlap between the two 
wave packets is very small and therefore very sensitive to the form of the tail of the 
wave function. At higher temperatures, however, the overlap is not small and the 
sensitivity to the form of the tail disappears. Consequently, (0) and (|53D become 
increasingly more accurate. 

At low temperatures, the second excited state is one of a pair of tunnel-split 
harmonic oscillator first excited states. The state with even parity does not contribute 
to the correlation function or to the approximate kink number formula because the 
matrix element of the position operator between the ground state and any even state 
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vanishes. The odd state, written in our double Gaussian approximation is, 

1 
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The expectation value of the energy in this state now follows from a straightforward 
calculation: 
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(60) 



In order to determine the values of 02 and ^2 one has in principle to carry out a 
constrained minimization of the energy by also enforcing the requirement that 
and \l/2 be orthogonal. From ( ^2]) and ([58|) this condition turns out to be 

2VL1VL2 - - 



exp 



nv2 



(61) 



^2 + 01 ^ \ Qi+Vt2 

For the range of temperatures we are interested in, the above condition is reasonably 
well satisfied if we take VLi = Q2 and (pi = 02, i-e., the left hand side of (61) is 
identically zero, while the right hand side is small (worst case of order 10^^). Therefore 
for this state we do not need to carry out a minimization of the energy as long as 
high accuracy is not required. In the low temperature limit. 
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The first term in (^) is just the energy of the first excited harmonic oscillator state 
as expected. The other terms represent anharmonic corrections and overlap contri- 
butions as usual. In Fig. V, the first three energies (as calculated in double Gaussian 
approximation) are plotted against the inverse temperature. While at intermediate 
temperatures A02 and Aqi are of similar magnitude, at low temperatures, A02 ^ Aqi 
which of course is just the WKB regime. It is clear that for values greater than /3 ~ 6, 
the WKB results should hold. 

The behavior of the energies of the ground and first two excited states (Fig. V) 
allows the identification of three qualitatively different regimes: (1) all the energies 
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lie above the classical barrier, (2) the ground state energy lies below the classical 
barrier height (at (3 = 1.734), and (3) the energy difference between the first two 
states becomes negligible in comparison with the energy difference between the ground 
and the second excited state (at /3 ~ 6). Our simulations confirm the theoretical 
expectations of the previous section that kinks cannot be identified in region (1), 
that there are kinks, but that the dilute gas approximation is invalid in region (2), 
and finally, that the dilute gas approximation is accurate in region (3). 

The classical single point field distribution function P[(f)] can be measured directly 
from our simulations. For the analogous quantum mechanical problem arising from 
the transfer operator method this is just the square of the ground state wave function 
\l/o. Results from the simulations are compared with "^^q at (3 = 2 and /5 = 4 in Fig. 
VI and are in reasonable agreement. The presence of kinks implies a double peak in 
P[4>] (the converse is false) while a single peak centered at the origin means that 
kinks and large amplitude thermal phonons can no longer be distinguished. From the 
simulations such a transition occurs at (3 ^ 1.7, in agreement with the theoretical 
calculation of when "^^q goes over from a double to single peaked distribution. (The 
double Gaussian method compares favorably with the numerical evaluation of the 
transfer operator in Ref. PD[, which predicts (3 ~ 1.8 as the transition point.) As 



expected, this is also the temperature (see Fig. 2) where the ground state energy 
crosses the classical barrier height. (A discussion of various methods to determine 
this characteristic temperature is given in Ref. ||24|| .) 



It is also apparent from Fig. VI that the double peaks in the distribution function 
move inward from the minimum of the classical potential as the temperature increases 
(eventually coalescing at (3 ~ 1.7). Physically this can be understood as nonlinear 
phonon corrections due to the fact that near each minimum, the potential is not 
symmetric under reflection around the minimum. 

The correlation function as determined by (|25| ) can of course be directly evaluated 
in the double Gaussian approximation. The correlation length is determined by the 
long distance behavior of the correlation function, and as is clear from (|38|) , one needs 
only the energy difference between the ground and first excited state to determine 
the correlation length. In Fig. Ill the correlation length as given by the double 
Gaussian and WKB approximations (see (|13|)) is compared with the results from the 
simulations. As expected, the double Gaussian prediction is borne out at intermediate 
temperatures (upto /3 ~ 4) whereas the WKB results are in agreement with the data 
at low temperatures (above /3 ~ 6). 

The kink number as determined by (^) requires the calculation of the matrix 
elements {0\(j)\n). The first two nonzero matrix elements as given by the double 



20 



Gaussian approximation are 

(O|0|l) = 

(O|0|2) = 




(63) 



where 
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As the energies for the first three states are aheady known we can now evaluate the 
kink densities as given by (^21). Theoretical estimates (for different choices of la and 
Ls) are compared with results from the simulation in Fig. II. The double Gaussian 
results for the kink density are entirely consistent with the numerical data aX (3 < 5.5 
(beyond this temperature the double Gaussian approximation underestimates the 
number of kinks due to the underestimation of Aqi). Due to the ambiguity in the 
concept of kink number discussed earlier, the significance of this comparison at high 
temperatures is not clear. 

The internal energy and the specific heat can be computed straightforwardly from 
the double Gaussian approximation. The behavior of the specific heat with temper- 
ature is of particular interest: it is known that at a temperature relatively close to, 
but somewhat below the characteristic temperature where kinks appear, the specific 
heat attains a maximum value. Though a real phase transition cannot occur in this 
system, it can be argued that this peak is a signal for the emergence of a new degree 
of freedom, in this case presumably, the kinks. While this peak has been seen in 
numerical calculations of the partition function it has not been unambiguously 
established that this feature is due to kinks. As remarked earlier, the ground state 
energy as given by the double Gaussian approximation breaks into two contributions: 
one given by an overlap contribution and arguably associated with kinks, the other 
associated with phonons. This feature allows one to analyze the nature of the peak 
in the specific heat. 

To begin, we recall the definitions. 
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and 2 

= -^^^ {PF) (69) 

for the internal energy U and the specific heat Cy. In our case, the free energy F 
has four contributions: one from the kinetic term in the original field Hamiltonian, 
one from the term given by the normalization of the functional integral, one from the 
overlap contribution, and one from the (potentially nonlinear) oscillations around the 
potential minimum of Hq. The first contribution stems from (^) giving rise to the 
free energy density 

(70) 

The other pieces come from (|52D. The contribution Fosc is just the first term of ( ^T] ) 
while Ftunn is the second term. The normalization contributes the Vq term. From 
([70|) and (|29|) we can immediately calculate the contribution to the internal energy 
density and the specific heat (per unit length) of the kinetic and measure terms. This 
gives just the results for a free field theory: 

^ i (71) 

C(») = 1. (72) 

Put another way, ([7l|) and (l7^) are the contributions from linear phonons. The 
nontrivial contributions are therefore included exclusively in Fosc and Ftunn- At low 
temperatures, the WKB result for the ground state energy may be used to show that 




2 ~ ^'^J 2 



1 „ 1 



e-'''". (73) 



vr 

In this regime Fosc does not contribute to the specific heat at all while the kink 
contribution is exponentially suppressed. The constant term due to linear phonons 
(fr^) is dominant. 

The situation can be dramatically different at intermediate temperatures. The 
specific heat as calculated from the double Gaussian approximation is plotted in Fig. 
VII: A prominent peak in the specific heat appears at /5 ~ 5.4. To understand whether 
this peak is due to kinks, the individual contributions from Fosc and F^nn are also 
plotted. There seems to be a delicate interplay between these two contributions. 
While the rise to the peak with increasing (3 is due to nonlinear phonons, just before 
the peak this contribution drops off steeply and eventually leads to a reduction in the 
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height of the peak. The overlap or kink contribution exhibits a slow initial decline as 
j3 increases, followed by a relatively sharp peak: it is this peak that is the dominant 
contribution at /5 ~ 5.4. We can therefore conclude that it is indeed the kinks that are 
responsible for the peak in the specific heat at this temperature. At larger values of 
/?, both contributions fall off smoothly to zero as expected from the low temperature 
result (|73|). 

We have also calculated the internal energy U . The nonlinear phonon and kink 
energies are plotted in Fig. VIII. There is a clear indication of an intricate energy 
sharing mechanism operating between the nonlinear phonons and kinks at intermedi- 
ate temperatures with kinks dominating the internal energy in the neighborhood of 
j3 ~ 5.4. Presumably this is due to the fact that at these temperatures kinks emerge 
as well defined localized objects and are created with relative ease. 

The results from the double Gaussian approximation as shown in Figs. VII and 
VIII are not quantitatively trustworthy for /5 > 5 as the approximation breaks down 
at that point. However, we expect them to be correct qualitatively: since this method 
underestimates the overlap at low temperatures, it in fact suppresses the kink con- 
tribution at lower temperatures. Therefore it is highly unlikely that the peak in the 
specific heat due to the kink contribution is an artifact of the approximation. A direct 
comparison of the theoretical results for the specific heat with numerical simulations 
is possible but requires high statistics. Work in this direction is in progress. 
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V. Discussion 

In summary, we have shown that the dilute gas/WKB approximation is excellent 
for j3 > 6 with no further phonon dressing of the bare kink energy beyond that already 
included in (0) and (|1^) at these low temperatures. In particular, we see no evidence 
of a reduction in the kink mass as suggested in some previous simulations At 
higher temperatures, the WKB analysis fails, though theoretical progress is possible 
with the double Gaussian technique. In this temperature regime, we have shown that 
the correlation length as measured in the simulations is in very good agreement with 
the theory and that the kink number density, while not being a well defined quantity, 
can still be understood at least qualitatively. The double Gaussian approximation 
also accurately predicts the onset of short-range order in the system as evidenced by 
the probability distribution function crossing over from a single to a double peaked 
distribution. 

The work reported here contains some of the largest simulations carried out to 
date on the model. The advantage of size is apparent: we have been able to go 
to low enough temperatures to unambiguously verify the WKB predictions for this 
model. Furthermore, the large system size enabled us to check for finite size effects 
and to be confident of their absence. 

The decomposition of the specific heat into two contributions via the double Gaus- 
sian approximation appears to demonstrate the existence of a nontrivial energy shar- 
ing interaction between kinks and nonlinear phonons in the neighborhood of the tem- 
perature where the specific heat is a maximum. At this temperature, {(3 ~ 5.4), kinks 
exist as well-defined objects (Fig. I) but the dilute gas/WKB theory is not valid. It 
would be interesting to see if an analytic study of the thermodynamics of kinks and 
phonons is possible in a phenomenological theory in this temperature regime. 

We expect to apply the double Gaussian variational method to other problems in 
the future, e.g., Sine-Gordon and $^ field theories in 1 + 1 dimensions. Work on the 
Langevin simulation method applied to calculate dynamical correlation functions is 
already in progress. 
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Figure Captions 

Fig. I: Sample field configurations, from top to bottom, at /3 = 2, /? = 4, /? = 5.5, 
and (3 = 8. Only a 1000 lattice unit sample of the total lattice size of 16384 is shown. 
Fig. II: Total number of kinks as a function of (3. Squares denote counts with a 
smoothed field [la = 8 lattice units) definition of kinks, diamonds for the zero-crossing 
counting method discussed in Sec. II, and the solid line starting at /5 = 4 is a fit to 
the WKB prediction (pH]). Also shown are three predictions (the three pairs of solid 
lines) for the kink number from the theoretical formula (|^) calculated in double 
Gaussian approximation. Three values for the averaging length were used, from top 
to bottom, /a=2,4, and 8 lattice units. For each value of la, there are two theoretical 
curves: (1) calculating ( ^2l) keeping only the ground and first excited states and (2) 
keeping the ground, and first and second excited states (the upper curves of each 
pair). The kink number is seen to strongly depend on la for (3 < Q. (The double 
Gaussian approximation breaks down at higher values of (3.) 

Fig. Ill: Field correlation length A as a function of (3. The double Gaussian (top 
curve) and WKB (bottom curve) predictions are compared with the numerical results. 
Squares are data points from the simulations. The crossover from the double Gaussian 
to the WKB range of validity happens around /3 ~ 5. 

Fig. IV: The Gaussian effective potential Vg (dashed lines) compared with Vdg 
(dot-dashed lines) at two temperatures {(3 = 2.5 and (3 = 4). The classical potential 
^ = -(1 - 07(2/32))07(2/52) from (0) is also shown (solid curves). For (3 < 3.12, 
Vg has a global minimum at the origin, while for (3 < 3.12 there are two degenerate 
global minima. For Vdg the origin is always a maximum and there are always two 
degenerate minima. The ground state energy calculated in the two approximations 
is the global minimum of Vg or of Vdg- The ground state energy is always lower in 
the double Gaussian approximation. 

Fig. V: The ground state, first excited, and the odd second excited state energies 
computed in double Gaussian approximation (the lowest, middle and uppermost solid 
lines) plotted against (3. The ground state energy has further been decomposed into 
the contributions from the "diagonal" (dot-dashed line) and "overlap" (dashed line) 
pieces. 

Fig. VI: The classical distribution function P[0] (solid lines) given by the simulation 
and the distribution \1/q^ from the double Gaussian approximation (dashed lines) 
plotted against for (3 = 2 and (3 = 4. The underestimation of the wavef unction 
overlap in the double Gaussian approximation is already visible at (3 = 4. 
Fig. VII: The specific heat Cy as calculated from the double Gaussian approximation 
plotted against (3 (top curve). The peak in the specific heat occurs at /5 ~ 5.4. The 
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contribution of the linear phonons to the specific heat is the constant value 1/A = 2. 
The individual contributions arising from the kinks (solid line) and nonlinear phonons 
(dot-dashed hne) are shown below. The WKB calculation for the kink contribution 
(dashed hne) is plotted from P — 6.5 onwards. 

Fig. VIII: The internal energy plotted against P in double Gaussian approximation 
(solid line). Kink (dashed line) and nonlinear phonon (dot-dashed line) contributions 
are shown separately (modulo an irrelevant constant energy shift). Note that kinks 
dominate over nonlinear phonons in the region where the specific heat has a maximum. 
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